m1 <- lm(I(log(Fst/(1-Fst)))~type,data= codingHCsnps,na.action=na.omit,subset=codingHCsnps$Fst>0&codingHCsnps$Fst<1)
tmp.x <-  m1$model
names(tmp.x) <- c("logit.fst","type")
m2 <- lm(logit.fst~type,data=tmp.x,subset=resid(m1)<10)

m1 <- lm(I(log(Fst/(1-Fst)))~type,data= trimmedHCsnps,na.action=na.omit,subset=trimmedHCsnps$Fst>0&trimmedHCsnps$Fst<1)
tmp.x <-  m1$model
names(tmp.x) <- c("logit.fst","type")
m2 <- lm(logit.fst~type,data=tmp.x,subset=resid(m1)<10)

mtmp <- trimmedHCsnps[trimmedHCsnps$Fst>0&trimmedHCsnps$Fst<1,]
mtmp <- mtmp[resid(m1)<10,]
mtmp$logit.fst <- log(mtmp$Fst/(1-mtmp$Fst))
library(lme4)

mm1 <- lmer(logit.fst~type+(1|cDNA),data=mtmp)

mtmp <- codingHCsnps[codingHCsnps$Fst>0&codingHCsnps$Fst<1,]
mtmp <- mtmp[resid(m1)<10,]
mtmp$logit.fst <- log(mtmp$Fst/(1-mtmp$Fst))

mm1 <- lmer(logit.fst~type+(1|cDNA),data=mtmp)

m1 <- lm(I(log(F/(1-F)))~type,data= trimmedHCsnps,na.action=na.omit,subset=trimmedHCsnps$F>0)


trimmedHCsnps$Fst <- apply(trimmedHCsnps,1,function(X){
	Fst.calc(p=as.numeric(X[grep('freq',names(X))])/100,sample.size=as.numeric(X[grep('cov',names(X))]))
	})

trimmedHCsnps$F <- apply(trimmedHCsnps,1,function(X){
	p <- as.numeric(X[grep('freq',names(X))])/100
	f.het <- (1-p)*p
	mean(f.het)
	})

###fsts etc
trimmedHCsnps$type <- factor("nonimmune",levels=c("nonimmune","immune"))
trimmedHCsnps$type[trimmedHCsnps$cDNA %in% tmp$id] <- "immune"
codingHCsnps <- trimmedHCsnps[trimmedHCsnps$cds,]
dnHCsnps <- trimmedHCsnps[trimmedHCsnps$dn,]

trimmedFst <- by(trimmedHCsnps[,c("Fst")], trimmedHCsnps$cDNA,mean)
trimmedF <- by(trimmedHCsnps[,c("F")], trimmedHCsnps$cDNA,mean)

tst <- by(trimmedHCsnps[,grep('freq',names(trimmedHCsnps))],trimmedHCsnps$cDNA,function(X){
	X <- X[!apply(X,1,function(X){any(is.na(X))}),]
	if(nrow(X)==0) return(rep(-1,ncol(X)))
	X <- X/100
	#pi <- rep(0,ncol(X))
	piX <- X*(1-X)
	colSums(piX)
	})

pi.df <- data.frame(cDNA = names(tst),pi.pop1=rep(0,length(tst)),pi.pop2=rep(0,length(tst)),
	pi.pop3=rep(0,length(tst)),stringsAsFactors=FALSE)
for(i in 1:length(tst)){
	pi.df[i,2:4] <- tst[[i]]
	}

tst2 <- by(codingHCsnps[,grep('freq',names(codingHCsnps))], codingHCsnps$cDNA,function(X){
	X <- X[!apply(X,1,function(X){any(is.na(X))}),]
	if(nrow(X)==0) return(rep(-1,ncol(X)))
	X <- X/100
	#pi <- rep(0,ncol(X))
	piX <- X*(1-X)
	colSums(piX)
	})

pi.df2 <- data.frame(cDNA = names(tst2),pi.pop1=rep(0,length(tst2)),pi.pop2=rep(0,length(tst2)),
	pi.pop3=rep(0,length(tst2)),stringsAsFactors=FALSE)
for(i in 1:length(tst2)){
	pi.df2[i,2:4] <- tst2[[i]]
	}
pi.df2$type <- factor("nonimmune",levels=c("nonimmune","immune"))
pi.df2$type[pi.df2$cDNA %in% imm.ids$id]  <- "immune"

tst3 <- by(dnHCsnps[,grep('freq',names(dnHCsnps))], dnHCsnps$cDNA,function(X){
	X <- X[!apply(X,1,function(X){any(is.na(X))}),]
	if(nrow(X)==0) return(rep(-1,ncol(X)))
	X <- X/100
	#pi <- rep(0,ncol(X))
	piX <- X*(1-X)
	colSums(piX)
	})

pi.df3 <- data.frame(cDNA = names(tst3),pi.pop1=rep(0,length(tst3)),pi.pop2=rep(0,length(tst3)),
	pi.pop3=rep(0,length(tst3)),stringsAsFactors=FALSE)
for(i in 1:length(tst3)){
	pi.df3[i,2:4] <- tst3[[i]]
	}
pi.df3$type <- factor("nonimmune",levels=c("nonimmune","immune"))
pi.df3$type[pi.df3$cDNA %in% imm.ids$id]  <- "immune"
